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[T , ! Abstract 

We present a systematic nonlinear theory of pattern selection for paramet- 
K^ ' ric surface waves (Faraday waves), not restricted to fluids of low viscosity. 

Q ' A standing wave amplitude equation is derived from the Navier-Stokes equa- 

(N' 

^^ ' tions that is of gradient form. The associated Lyapunov function is calculated 

l> 

5^ ' for different regular patterns to determine the selected pattern near thresh- 

ij^ • old. For fluids of large viscosity, the selected wave pattern consists of parallel 






C^ 



stripes. At lower viscosity, patterns of square symmetry are obtained in the 
capillary regime (large frequencies). At lower frequencies (the mixed gravity- 



/\ • capillary regime), a sequence of six-fold (hexagonal), eight-fold, . . . patterns 



are predicted. The regions of stability of the various patterns are in quan- 
titative agreement with recent experiments conducted in large aspect ratio 
systems. 
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Parametrically driven surface waves (also known as Faraday waves) appear on the free 
surface of fluid layer which is periodically vibrated in the direction normal to the surface 
at rest. Above a certain critical value of the driving amplitude, the planar surface becomes 
unstable to a pattern of standing waves |1[] . If the viscosity of the fluid is large, the bifurcating 
wave pattern consists of parallel stripes. At lower viscosity, patterns of square symmetry 
are observed in the capillary regime (large frequencies) |^ . At lower frequencies (the mixed 
gravity-capillary regime), hexagonal, eight-fold, and ten-fold patterns have been observed 
P-§[|. We present a weakly nonlinear analysis of the equations governing fluid motion that 
predicts stationary wave patterns with these symmetries. Their boundaries of stability agree 
quantitatively with experiments. [^Q 

Pattern selection in confined geometries can be often understood in terms of the spatial 
modes of the base state that become linearly unstable and the geometry of the system. 
Extended systems, on the other hand, allow a richer nonlinear competition of linearly un- 
stable modes, in part due to the restoration of some of the symmetries of the original system 
that had been broken by the boundaries. The distinction can be further quantified by in- 
troducing the coherence length of the pattern ^, the linear size of the system L, and the 
characteristic wavelength of the pattern I/Zcq- Directional solidification from the melt is a 
typical example of a system in the limit ^o ^ l/^o ^ L, fact that follows from the ex- 
tremely flat neutral stability curve 0. Experiments often show a narrow distribution of 
wavelengths in the stationary state, the prediction of which has proven elusive. Attempts 
at deriving amplitude equations valid near onset have not been successful precisely because 
of the condition ^q ^ l/^o- A typical and widely studied example of the intermediate range 
concerns Rayleigh-Benard convection in large aspect ratio cells I/Zcq i4 ^ '^ L. Although 
amplitude equations predict the existence of stable parallel rolls above onset, such a state 
is not generically observed in fluids of low or moderate Prandtl number except under care- 
fully prepared initial conditions. Instead, a spatio-temporally chaotic state emerges that has 
been termed spiral defect chaos [0. The role that long wavelength modes (or mean flow) 
play in the development of such a state is still a matter of research. Faraday waves, on the 



other hand, are a prototypical case of a system with a large coherence length ^ ^ l/A;o, and 
therefore one would expect that amplitude equations would allow quantitative predictions 
of pattern selection near onset. Furthermore, the physical system and the experimental con- 
ditions are completely determined by a few independent parameters that can be obtained 
with reasonable accuracy. 

Progress in deriving suitable amplitude equations for Faraday waves has proved difficult. 
Since Faraday waves are almost a Hamiltonian system (weakly dissipative) , most analyses 
are based on a hamiltonian description for the ideal (inviscid) limit, and treat viscous or 
dissipative effects as a perturbation [§,0. The central question that arises in this case 
concerns the origin of any nonlinear saturation of the wave pattern. First, general symmetry 
considerations for Hamiltonian systems ( "chiral symmetry" ) prohibit cubic terms in standing 
wave amplitude equations [||,|l3]- In addition, it has been argued that linear viscous terms in 



the original equations governing Faraday waves only contribute linear damping terms to the 
amplitude equations ||^. These observations have contributed to the belief that saturation in 
near-Hamiltonian systems generically occurs through either nonlinear (cubic) viscous terms 
or fifth order hamiltonian terms, both of which are very difficult to obtain explicitly [Q. 
The issue, however, has remained somewhat controversial for the case of Faraday waves 
IP, and we mention, for example, two recent studies that considered non-dissipative cubic 



terms in the amplitude equation. Miiller |11| has shown that the generic form of a cubic 



amplitude equation would allow patterns of standing waves with the observed symmetries, 
given suitable choices of the coefficients. A more recent analysis |1^ explicitly addressed 



these issues in the limit of small fluid viscosity. A set of quasi-potential equations was 
derived by considering the rotational flow within a small viscous layer near the free surface, 
and assuming potential flow in the bulk. A standing wave amplitude equation was obtained 
with cubic order terms given explicitly. 

Although this latter model predicted bifurcations to square and higher-symmetry pat- 
terns, it did rely on an uncontrolled approximation concerning nonlinear viscous terms. As 
a consequence, its region of validity is difficult to asses. In particular, it failed to yield the 



observed stripe patterns at intermediate and large viscosities. We present here a general 
calculation that overcomes these difficulties, and that leads to the experimentally observed 
regular (periodic or quasi-periodic) standing wave patterns above onset. As part of the 
derivation, we also obtain an analytical expression for the linear threshold of instability, 
which was previously known only numerically |T^. Not being confined to small viscosities 



is also important for comparison with experiments since in this case it is easier to achieve 
the large aspect ratio limit, and hence to study pattern formation without the influence of 
side walls. 

We consider an incompressible viscous fluid under vertical vibration /cos (cut) of ampli- 
tude /, and angular frequency uj. The fluid at rest has a free surface at z = 0, extends to 
z = — oo, and is unbounded in the x-y direction. The equation governing fluid motion is 

dtu + (u • V)u = — Vp + uV^u + G(t)e,, 
P 

with u the velocity field, p the pressure, p and z/ the density and kinematic viscosity of the 

fluid respectively, and G(t) = —g — |/(e*'^* + e"""^*) the effective gravity. For / below the 

threshold of instability, the base state is u = and p = pG{t)z. We ffist eliminate the 

explicit dependence on the pressure by taking —V x Vx to obtain 

dtV^u - z/V^ V^u = V X V X (u ■ V)u. 

Here the continuity equation V ■ u = has also been used. 

The position of the free surface is denoted by ({x,y), the unit normal n = 
{—dx(, —dy(, 1), and the two tangential unit vectors are ti = (1, 0, d^C) and t2 = (0, 1, dy(). 
Besides the null conditions at z = — oo, there are three boundary conditions to be satisfied 
at the free surface, 

dtC + (u ■ Vh)C = uj\z=<: 

t^ ■ T ■ n|^=^ = 0, m = l,2 



with Vh = ^xdx + ^ydy, T the stress tensor with components T^- = [—p — pG(t)z]Sij + 
pv{djUi + diUj), a the surface tension, and 2H the mean curvature of the free surface |T^. 
First we consider the hnear stabihty of a subharmonic standing wave, 

Wo = cos{kx) Y. e^''^'/^wi{z)Aj + c.c, 

where Wq is the ^-component of the velocity field, and a similar expansion for (q. Substitution 
into the linearized equation of motion, (9fV^ — z/V^V^) wq = 0, and the linearized kinematic 
and tangential stress boundary conditions, dtCo — wq = and (V|^ — d1)wQ = 0, we find 

wi{z) = vi^ + q])e^'' - 2uk^e''^\ 

with q^- = k"^ + jiuj jlv. Note that the boundary conditions at 2; = C have been expanded 
around 2; = 0. (Generalization to a finite fluid depth is straightforward: the term d\e~^''^ + 
(^2e~*^ will also be in Wo(-z) with d\ and ^2 to be determined by the additional boundary 
conditions at the bottom.) 

The critical amplitude /o is determined by the linearized normal stress boundary condi- 
tion, which is, after using the momentum equation to eliminate po^ 

2pv\l\ d^ Wq - pdtd^wo + pvV'^d^WQ + pgV]jCo 

+ \pf (e*-* + e-*-*) ViCo - aV^V^Co = 0. 

By substituting wq and Co into the above equation, we find for each harmonic e-'*'^*/^, 

iJiAi - fAl - /A3 = 
H,A, - fA, - /A = 
H,A, - fAs -fAj = 

with Hj = Ipu'^ Aqjk^ — k{q^ + k'^Y — pgk"^ — ak^\ /\pk'^. By truncating the set of equa- 
tions dll) at some y4„, the system can be solved numerically as an eigenvalue problem. This 
is indeed what was done by Kumar and Tuckerman ||13|. However we observe that after 



truncation at An, v4„ = fAn-2/Hn, An-2 = /^n-4/(-f^n-2 — jj-), ■ ■ ■■ Therefore the set of 
equations can be rewritten as 

Hi - ^ ^\f, J A, - fAl ^ H,ik,f)A, - fAl = 0, 

so that for a given wavenumber k, the threshold of instabihty /q is given by 

fo = \H,{k, fo)\. 

The complex amplitude Aj can be recursively obtained from Eq. (|l]) up to a real factor. For 
an infinite system the critical wavenumber /conset is the wavenumber that corresponds to the 
lowest value of /q. 

Consider first the limit of low viscosity, and define Uq = uj/2 as the time scale, and ko 
from cUq = gkQ + crkQ/p as the length scale, and the dimensionless variables k = k/kQ,'y = 
2i'kQ/ujQ, G = gko/uQ, S = akH puj^, and A = foko/AuQ. For a damping coefficient 7 <C 1 
and k near fconset, Aonset can be given explicitly as, 

A 1 3/2 5/2 

Aonset = 7-27 +8(3_2G')^ +■■■' 
with < G < 1 by definition. The dimensionless critical amplitude is proportional to 
7 — ^7'^''^ at small 7. While previous low damping calculations |]8|,p!2[| only used the linear 
term to determine the location of the threshold, the first correction — 17^'^ can be a sizable 
contribution (e.g., a 15% difference at 7 = 0.1). We also note that a similar calculation for 
the damped Mathieu equation leads to a threshold 7 + 37^/64 + 0(7^), in which the first 
correction term is of a different order and has a different sign. 

To derive the amplitude equation we use the multiple scale approach |TB|. The solvability 



condition in this case arises from the boundary conditions, not from the equation of motion 
as in most other cases. The velocity field is expanded as, 

u = e^/^uo + eui + e^/^Ua + • • • , 

with e = (/ — /o)//o, and similarly for p and (. Near threshold, i.e., for e <^ 1, we 
separate fast and slow time scales: T = et; dt ^ dt + edr- Spatial slow scales are not 



included because only regular patterns are considered here. At order e^/^ we recover the 
linear solution discussed above. Because we are interested in standing wave patterns with 
different symmetries, the solution at this order is written as a linear combination of waves 
with wavevectors k^ of magnitude fconset in different directions on the x-y plane, 

«7o = E cos(k^ ■ v)B^{T) J2 e^"^'/^wiiz)e, + c.c. 

m j=l,3,5,-- 

Here Bm{T) are the real wave amplitudes, functions only of the slow time scale T, and Cj is 
Aj found in Eq. (p. 

At order e the equation of motion for wi becomes 

(dtV^ - z/V^V^) wi = [V X V X (uo ■ V)uo],. (2) 

The solution Wi contains terms of the form cos ((k^ ± k„) ■ r) that incorporate couplings 
through stable modes, and that will contribute to the coefficients of the cubic terms later 
in the expansion. The particular solution wip is obtained by integrating Eq. (H), with the 
homogeneous solution wi^ chosen so that the boundary conditions are satisfied. Both the 
equation of motion and boundary conditions at this order become very complicated. In 
order to find Wi in practice, we have developed a symbolic manipulation program specific 
to this case, and found the solution on a computer. 
At order e^^^ the equation of motion becomes 

{dtV' - z/VV^) W2 (3) 

= -drV^wo + {V X V X [(uq ■ V)ui + (ui • V)uo]}, . 

Only terms proportional to cos(ki ■ r) need to be considered in the solution for W2 and (2, 
W2 = cos(ki ■ r) E e^*""*^^ [^j + {cije^' + h^^'"'") Cj 

i=l,3,5,- 

C2 = cos(ki-r) E e^'^^'/'C,. 

J=l,3,5,- 

Here Ej{z) comes from the direct integration of Eq. (^), and ajC^^+bjC'^^^ is the homogeneous 
solution that has the same form as the linear solution. 



Using the kinematic and tangential stress boundary conditions at order e"!"^ we find (again 
with the symbohc manipulation program) a^ and hj. The solution w^ is finally inserted into 
the normal stress boundary condition at order e"!"^ to yield a system of equations for Cj 
which has the same left-hand side as Eq. (|l]) but with nonzero right hand side. Solving for 
Cj just like in Eq. (|l|), we obtain 

HiC^ - Wl = F, 

with F a function of the amplitude B^- Since /o = \Hi\, by requiring a nontrivial solution 
for Ci we obtain the solvability condition FH^ + -F*/o = 0, which yields a standing wave 
amplitude equation, 

JJD 

-^ = aB,- goBf - ^ g{emi)BlB,, (4) 

with 6ml the angle between k^, and ki, and explicit expressions for the coefficients. Equation 
(I) is of gradient form, and can be derived from a Lyapunov function 

^ = -lc.EBl + \Y.J:9{0n.n)BlBl (5) 

rri m n 

which may be used to find the preferred pattern near threshold |I]]. For regular patterns of 
A^ standing waves, k^ form a regular polygon, and the Bm are constant. 

We next turn to a comparison between the selected patterns predicted by Eq. (|), and 
two recent sets of systematic experimental surveys involving large aspect ratio systems, both 
of which aiming at addressing the issue of pattern selection over a wide range of experimental 
parameters @J^. Binks and van de Water have focused on a low viscosity fiuid [|1^, a 



large aspect ratio cell, and a layer depth much larger than the wavelength. When the driving 
frequency is decreased from 45Hz, a transition from a A^ = 2 square pattern to a A^ = 3 
hexagonal pattern was observed at approximately 35Hz, and to a quasi-periodic A^ = 4 
eight-fold pattern at approximately 29Hz [|I^. Our prediction for these transitions based on 
Eq. (^ are 35.4 Hz and 28.7 Hz respectively. These results are also in good agreement with 
the earlier weak damping calculation [Q that predicted the same transitions at frequencies 
of 32.8 Hz and 27.9 Hz respectively. 

8 



A large aspect ratio experiment involving fluids of various viscosities has been carried 
out by Kudrolli and Gollub [^. Although the fluid depth (0.3 cm) is smaller than the 
wavelength in the experiment (1-3 cm), the comparison is still illuminating. Figure |1| shows 
the symmetry of the preferred patterns predicted by our calculations in the parameter space 
defined by the viscosity of the fluid and the driving frequency (with p = 0.95g/cm^ and 
a = 20.6dyne/cm), and the experimentally observed patterns. Stripe patterns are preferred 
at high viscosity, whereas at low viscosity, hexagons (at low frequency) and squares (at high 
frequency) are observed. We also show the small region in which a sequence of quasiperiodic 
patterns are expected to be selected. The experimental results by Kudrolli and Gollub 
are shown as the symbols in the figure. We note an excellent agreement in the regions in 
which stripes, squares, and hexagons are observed, specially in view of the shallow fluid 
depth involved in the experiment. The shallowness of the layer probably accounts for the 
observation of a hexagonal pattern at z/ = Icm^/s and low frequency, and not observing a 
quasiperiodic pattern for u = 0.04cm^/s and / = 27Hz. As noted above, the experiments 
by Binks and van de Water [^ did probe this latter region in a deep fluid layer, with their 
results agreeing with our predictions. 

In summary, we have presented a nonlinear theory for Faraday waves in viscous fluids 
with no assumptions or approximations other than those inherent to the multiscale expan- 
sion. A set of standing wave amplitude equations has been obtained that is of gradient 
form. Minimization of the associated Lyapunov function leads to determination of the pre- 
ferred pattern near threshold. The predicted patterns are in excellent agreement with recent 
experiments in large aspect ratio systems involving a range of fluid viscosities and driving 
frequencies. 
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FIG. 1. Preferred patterns in viscosity-driving frequency space. Symbols represent the ex- 
perimental results. x=stripe, □=square, and A=hexagon. Alternating x and □ indicate 
mixed-stripe-square states. 
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